

In [1]:

    
name = "Shehzad Qureshi"
jhed_id = "squresh6"
if name == "Student Name" or jhed_id == "sname1":
    raise Exception( "Change the name and/or JHED ID...preferrably to yours.")



In [2]:

    
from IPython.core.display import *
from Queue import PriorityQueue
import math

Search is often used for path-finding in video games. Even though the characters in a video game often move in continuous spaces, it is trivial to layout a "waypoint" system as a kind of navigation grid over the continuous space. Then if the character needs to get from Point A to Point B, it does a line of sight (LOS) scan to find the nearest waypoint (let's call it Waypoint A) and finds the nearest, LOS waypoint to Point B (let's call it Waypoint B). The agent then does a A* search for Waypoint B from Waypoint A to find the shortest path. The entire path is thus Point A to Waypoint A to Waypoint B to Point B.

We're going to simplify the problem by working in a grid world. The symbols that form the grid have a special meaning as they specify the type of the terrain and the cost to enter a grid cell with that type of terrain:

token terrain cost . plains 1

forest 3 ^ hills 5 ~ swamp 7 x mountains impassible </code>

When you go from a plains node to a forest node it costs 3. When you go from a forest node to a plains node, it costs 1. You can think of the grid as a big graph. Each grid cell (terrain symbol) is a node and there are edges to the north, south, east and west (except at the edges).

For this problem, the agent will start in (0, 0) or the upper left hand cell and is looking for the lowest cost path to the goal at (26,26) which is in the lower right. Of course, your code should work for any start and goal coordinates on the map. You need to implement A* search to get the agent from start to goal.

The code below will read the world file. Please always put resource files in the same directory/folder as your notebook and use relative references (I don't have your C drive on my Mac!).



In [3]:

    
with open('world.txt', 'r') as f:
    map_file = [x for x in f.readlines()]
f.closed









    Out[3]:





True

First, we're going to convert the "world" into a map represented as a list of lists of character strings.



In [4]:

    
world = []
for line in map_file:
    line = line.strip()
    if line == "": continue
    world.append([x for x in line])

The map is now a list of lists. Keeping with regular x (across) and y (down) coordinates, you now have to access a point by first referencing y and then referencing x. So that world[2][3] is (3, 2).

We're going to need a few helper data structures such as costs as a dict(ionary):



In [5]:

    
costs = { '.': 1, '*': 3, '^': 5, '~': 7}

and a list of offsets for neighbors which is really a movement model:



In [6]:

    
reachable = [(0,-1), (1,0), (0,1), (-1,0)]

A* Search

This algorithm was implemented using the psuedocode in the textbook as a guide. A* search requires the use of a Priority Queue, however, python has no data structure that can serve as a priority queue and do efficient membership testing as well. In the end, helper functions were implemented to assist with the "queue" implementation, which really is just a list. Using a list isn't the most efficient way of doing this but it certainly is the easiest.

The heurisitic or evaluation function used to supplement the path cost was a straight-line distance from the node to the goal. The higher the value, the least preferable the action. The straight-line distance was calculated using the distance formula, i.e. the distance between two points.

function frontier_exists

Frontier helper function that checks to see if a position/coordinate exists in the frontier. If it does, the (cost, position) tuple is returned otherwise None is returned.



In [7]:

    
def frontier_exists(frontier, position):
    exists = any(x for x in frontier if x[1] == position)
    if exists is False:
        return None
    else:
        return [x for x in frontier if x[1] == position][0]

function frontier_empty

Frontier helper function that checks if the frontier is empty. Returns True or False.



In [8]:

    
def frontier_empty(frontier):
    return len(frontier) == 0

function frontier_insert

Frontier helper function that inserts a new (cost, position) tuple into the frontier. If the position already exists an exception is raised.



In [9]:

    
def frontier_insert(frontier, position, cost):
    if frontier_exists(frontier, position) is None:
        frontier.append((cost, position))
    else:
        raise ValueError("Position already exists!")

function frontier_update

Frontier helper function that updates a (cost, position) tuple already in the frontier. If the position does not exist an exception is raised. Iteration is required instead of using the other helper functions because python does not allow modifying tuples. This, we need the index of the specified position to overwrite it with a new (cost, position) tuple.



In [10]:

    
def frontier_update(frontier, position, cost):
    for index, item in enumerate(frontier):
        if item[1] == position:
            frontier[index] = (cost, position)
            return
    raise ValueError("Position does not exist!")

function frontier_pop

Frontier helper function that removes the first item in the frontier which will have the lowest cost value. The actual list itself is not sorted.



In [11]:

    
def frontier_pop(frontier):
    item = sorted(frontier)[0]
    frontier.remove(item)
    return item

function reconstruct_path

Algorithm helper function that reconstructs the proper path taken given all the nodes explored.

In order to properly reconstruct the path, we need to make certain assumptions:

Movement only occurs from one position to its adjacent position, i.e. similar to Hamming distance
Movement to its adjacent position gives a path cost of 1.

The path explored is iterated over in reverse order (i.e. from goal to start). If the distance between the current position and the previous position does not give a cost of 1, we can assume that the search algorithm branched to a different path to try and find a lower cost alternative.

The output of this function will yield a path explored from start to goal.



In [12]:

    
def reconstruct_path(path):
    new_path = []
    current_node = path.pop()
    new_path.append(current_node)
    while len(path) > 0:
        parent_node = path.pop()
        distance = get_distance_cost(current_node, parent_node)
        if distance == 1:
            new_path.append(parent_node)
            current_node = parent_node
    return [x for x in reversed(new_path)]

function get_distance_cost

Algorithm helper function that calculates the straight-line distance between two points.



In [13]:

    
def get_distance_cost(p1, p2):
    return math.sqrt((p2[0] - p1[0])**2 + (p2[1] - p1[1])**2)

function get_path_cost

Algorithm helper function that retrieves the path cost of reaching that position based on the the world map and cost value provided by the user.



In [14]:

    
def get_path_cost(world, costs, position):
    x = position[0]
    y = position[1]
    cost = costs[world[y][x]]
    return cost

function get_total_cost

Algorithm helper function that gets the total cost incurred for reaching that position. This includes the cost to get to this position and the hypothetical cost to get to the goal.



In [15]:

    
def get_total_cost(world, costs, position, goal):
    distance = get_distance_cost(position, goal)
    cost = get_path_cost(world, costs, position)
    return cost + distance

function possible_actions

Algorithm helper function that determines the potential possible actions from the current/given position. We need to make sure that the next potential is:

Within the bounds of the map
Has an associated cost function

A list of potential actions are returned. Note that we do not check if the potential action is in the explored or frontier list; that is done in the search algorithm function.



In [16]:

    
def possible_actions(world, costs, reachable, position):
    actions = list()
    for dx, dy in reachable:
        x = position[0] + dx
        y = position[1] + dy
        if not is_legal_position(world, costs, (x, y)): continue
        actions.append((x, y))
    return actions

function is_legal_position

Algorithm helper function that determines if a position is legal. A position is legal if it is within the bounds of the map and it has a cost associated with it. In our scenario, an unpassable region denoted by 'x' would not be legal.



In [17]:

    
def is_legal_position(world, costs, position):
    bound_x = len(world[0]) - 1
    bound_y = len(world) - 1
    x, y = position
    if x < 0 or x > bound_x: return False
    if y < 0 or y > bound_y: return False
    if world[y][x] not in costs: return False
    return True

function a_star_search

This is the main entry point for the A* search algorithm. It takes in a world, costs parameters, movement model and a start and goal coordinates within the world. A* search iterates through the world, starting at the goal, using a path cost and a heuristic model previously defined to determine the lowest cost path to the goal. This code was adapted from the psuedocode of Uniform Search Cost on page 86 of Artificial Intelligence: A Modern Approach Third Edition by Russell et al.

There wasn't a need to keep a path_cost variable since the frontier contained updated path costs anyway. A list was used to keep track of the order of positions explored, since the explored data structure does not maintain order. The algorithm is otherwise similar to the psuedocode presented in the book.



In [18]:

    
def a_star_search(world, costs, reachable, start, goal):
    if not is_legal_position(world, costs, start) or not is_legal_position(world, costs, goal):
        print("Start or Goal is illegal")
        return []
    frontier = [(0, start)]
    explored = set()
    path = []
    while True:
        if frontier_empty(frontier) is True:
            print("Could not find goal")
            return []
        current_cost, node = frontier_pop(frontier)
        path.append(node)
        explored.add(node)
        if node == goal:
            return reconstruct_path(path)
        for action in possible_actions(world, costs, reachable, node):
            item = frontier_exists(frontier, action)
            cost = get_total_cost(world, costs, action, goal) + current_cost
            if action not in explored and item is None:
                frontier_insert(frontier, action, cost)
            elif item is not None:
                if cost < item[0]:
                    frontier_update(frontier, action, cost)

Test

Finally, our test function which performs an A* Search from the starting point (0, 0) to the goal (26, 26). It outputs a list of coordinates/positions of the path taken from the start to the goal.



In [19]:

    
a_star_search(world, costs, reachable, (0, 0), (26, 26))









    Out[19]:





[(0, 0),
 (0, 1),
 (0, 2),
 (0, 3),
 (0, 4),
 (0, 5),
 (0, 6),
 (0, 7),
 (0, 8),
 (0, 9),
 (0, 10),
 (0, 11),
 (0, 12),
 (1, 12),
 (2, 12),
 (3, 12),
 (4, 12),
 (4, 13),
 (5, 13),
 (5, 14),
 (6, 14),
 (6, 15),
 (6, 16),
 (6, 17),
 (7, 17),
 (8, 17),
 (9, 17),
 (10, 17),
 (11, 17),
 (12, 17),
 (13, 17),
 (14, 17),
 (15, 17),
 (15, 18),
 (16, 18),
 (17, 18),
 (18, 18),
 (19, 18),
 (20, 18),
 (21, 18),
 (22, 18),
 (23, 18),
 (23, 19),
 (23, 20),
 (23, 21),
 (23, 22),
 (23, 23),
 (23, 24),
 (23, 25),
 (24, 25),
 (24, 26),
 (25, 26),
 (26, 26)]

Summary

A* is an interesting algorithm which I have repeatedly heard about but am finally learning about it. I was surprised to learn that the algorithm is actually similar to a Breadth First Search and not a Depth First Search as I previously assumed. A BFS algorithm is usually pretty exhaustive which is what I suppose makes this algorithm complete and optimal, although I was shocked to learn that a DFS is neither complete nor optimal. I will have to investigate why later, time-permitting.

I initially implemented the Frontier data structure as a class, which really just contained all the helper functions within itself and operated on a list. I was able to refactor it out of the OOP style into an iterative style although I must admit that I may have liked the class implementation better. Purely for aesthetic reasons though.

I did spend a lot of time muddling my way through IPython, which the majority of it on debugging. It seems that there really isn't a good way of debugging in the notebook itself (none that I can find) without resorting to scattering print statements everywhere. I ended up using a regular python shell and copy/pasted the code in so that I can step through the program and find the problem quicker.

I did want to investigate the ordering of the movement model. Specifically, what would happen if I had a "right" movement before the "down" movement? Would this change the results of the algorithm? The function possible_actions iteratively checks the potential actions dependent on how the movement model is defined.



In [20]:

    
reachable = [(0,-1), (1,0), (0,1), (-1,0)] # original: up, right, down, left
test_reachable_order_1 = a_star_search(world, costs, reachable, (0, 0), (26, 26))
reachable = [(0,1), (1,0), (0,-1), (-1,0)] # down, right, up, left
test_reachable_order_2 = a_star_search(world, costs, reachable, (0, 0), (26, 26))
test_reachable_order_1 == test_reachable_order_2









    Out[20]:





True



In [21]:

    
reachable = [(1,0), (0,1), (-1,0), (0,-1)] # right, down, left, up
test_reachable_order_3 = a_star_search(world, costs, reachable, (0, 0), (26, 26))
test_reachable_order_1 == test_reachable_order_3









    Out[21]:





True

So the order in which we specify the movement model really does not matter! I suppose that scores one for optimality?

I did want to actually see my data visually. The path should be plottable assuming I reverse the y-axis and plot it properly. This helper function helps to setup the basic plot.



In [22]:

    
def setup_plot():
    fig, ax = plt.subplots()
    ax.set_xlim((-1, len(world[0])))
    ax.set_ylim((-1, len(world)))
    plt.gca().invert_yaxis()
    return fig, ax

Now we setup the map background. '~' does not exist in matplotlib so we will use a diamond instead. 'x' has a higher alpha value to stand out more.



In [23]:

    
def setup_map(ax, world):
    c = dict()
    c['.'] = ('.', 0.2)
    c['x'] = ('x', 0.7)
    c['^'] = ('^', 0.2)
    c['*'] = ('*', 0.2)
    c['~'] = ('d', 0.2)
    for y, _y in enumerate(world):
        for x, _x in enumerate(_y):
            params = c[world[y][x]]
            ax.plot(x, y, params[0], color='black', alpha=params[1])

Another helper function to calculate the total path cost.



In [24]:

    
def get_total_path_cost(world, costs, path):
    tcost = 0
    for node in path:
        x, y = node
        tcost = tcost + costs[world[y][x]]
    return tcost

Some imports to help us.



In [25]:

    
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

Re-run the path and plot the result!



In [26]:

    
path = a_star_search(world, costs, reachable, (0, 0), (26, 26))
print("Total path cost: {0}".format(get_total_path_cost(world, costs, path)))
x, y = zip(*path)

fig, ax = setup_plot()
fig.set_size_inches(10, 10)
setup_map(ax, world)
ax.plot(x, y)
plt.show()









    



Total path cost: 169

So it does look like the path is calculated properly. It does not go through any impassable regions but I do wonder if this indeed is the optimal path through the world. It is entirely possible to move down from (0, 0) through (0, 22) and then make your way across, but why does it not do this? I will assume that the heuristic (straight-line distance) will start to favor crossing the map in the middle of the map, which is what I see. What happens if we don't use the heuristic?



In [27]:

    
def best_first_search(world, costs, reachable, start, goal):
    if not is_legal_position(world, costs, start) or not is_legal_position(world, costs, goal):
        print("Start or Goal is illegal")
        return []
    frontier = [(0, start)]
    explored = set()
    path = []
    while True:
        if frontier_empty(frontier) is True:
            print("Could not find goal")
            return []
        current_cost, node = frontier_pop(frontier)
        path.append(node)
        explored.add(node)
        if node == goal:
            return reconstruct_path(path)
        for action in possible_actions(world, costs, reachable, node):
            item = frontier_exists(frontier, action)
            cost = get_path_cost(world, costs, action) + current_cost
            if action not in explored and item is None:
                frontier_insert(frontier, action, cost)
            elif item is not None:
                if cost < item[0]:
                    frontier_update(frontier, action, cost)

Re-plot.



In [28]:

    
path = best_first_search(world, costs, reachable, (0, 0), (26, 26))
print("Total path cost: {0}".format(get_total_path_cost(world, costs, path)))
x, y = zip(*path)

fig, ax = setup_plot()
fig.set_size_inches(10, 10)
setup_map(ax, world)
ax.plot(x, y)
plt.show()









    



Total path cost: 109

Looks much better! So the problem is with my heuristic! Let's try again using only heuristics, or greedy search.



In [29]:

    
def greedy_search(world, costs, reachable, start, goal):
    if not is_legal_position(world, costs, start) or not is_legal_position(world, costs, goal):
        print("Start or Goal is illegal")
        return []
    frontier = [(0, start)]
    explored = set()
    path = []
    while True:
        if frontier_empty(frontier) is True:
            print("Could not find goal")
            return []
        current_cost, node = frontier_pop(frontier)
        path.append(node)
        explored.add(node)
        if node == goal:
            return reconstruct_path(path)
        for action in possible_actions(world, costs, reachable, node):
            item = frontier_exists(frontier, action)
            cost = get_distance_cost(action, goal) + current_cost
            if action not in explored and item is None:
                frontier_insert(frontier, action, cost)
            elif item is not None:
                if cost < item[0]:
                    frontier_update(frontier, action, cost)



In [30]:

    
path = greedy_search(world, costs, reachable, (0, 0), (26, 26))
print("Total path cost: {0}".format(get_total_path_cost(world, costs, path)))
x, y = zip(*path)

fig, ax = setup_plot()
fig.set_size_inches(10, 10)
setup_map(ax, world)
ax.plot(x, y)
plt.show()









    



Total path cost: 155

We still run into the same problem, i.e. the algorithm favors cutting across the map at the middle because the distance to the goal is smaller at that point. So what happens if we restrict the straight-line distance to a certain range?



In [31]:

    
def get_distance_cost(p1, p2):
    dist = math.sqrt((p2[0] - p1[0])**2 + (p2[1] - p1[1])**2)
    if dist <= 3:
        return dist
    else:
        return 3



In [32]:

    
path = a_star_search(world, costs, reachable, (0, 0), (26, 26))
print("Total path cost: {0}".format(get_total_path_cost(world, costs, path)))
x, y = zip(*path)

fig, ax = setup_plot()
fig.set_size_inches(10, 10)
setup_map(ax, world)
ax.plot(x, y)
plt.show()









    



Total path cost: 123

Much better, although I suspect that by restricting the heuristic we're moving further away from A* and closer to best-first.

So what if we took the movement model and added diagonals?



In [33]:

    
reachable = [(0,-1), (1,0), (0,1), (-1,0), (1,1), (-1,1), (-1,-1), (1,-1)]

Reset distance calculator.



In [34]:

    
def get_distance_cost(p1, p2):
    return math.sqrt((p2[0] - p1[0])**2 + (p2[1] - p1[1])**2)

Our path reconstruction needs to be modified because diagonals do not have unit cost.



In [35]:

    
def reconstruct_path(path):
    new_path = []
    current_node = path.pop()
    new_path.append(current_node)
    while len(path) > 0:
        parent_node = path.pop()
        distance = get_distance_cost(current_node, parent_node)
        if distance == math.sqrt(2) or distance == 1:
            new_path.append(parent_node)
            current_node = parent_node
    return [x for x in reversed(new_path)]

Re-run the search again with diagonal movements this time.



In [36]:

    
path = a_star_search(world, costs, reachable, (0, 0), (26, 26))
print("Total path cost: {0}".format(get_total_path_cost(world, costs, path)))
x, y = zip(*path)

fig, ax = setup_plot()
fig.set_size_inches(10, 10)
setup_map(ax, world)
ax.plot(x, y)
plt.show()









    



Total path cost: 140

So the path looks better albeit a bit more jagged than I would expect. I suspect that the path reconstruction is at fault because it would need to predict whether to pick a diagonal next or a lateral movement. This is not as easy as inspecting the next element since it might decide to move lateral or diagonal twice, thereby throwing off the reconstruction. Time-permitting, I would like to investigate better path reconstruction. And just for fun, we can try best first search and greedy search again to see if they are any better.



In [37]:

    
path = best_first_search(world, costs, reachable, (0, 0), (26, 26))
print("Total path cost: {0}".format(get_total_path_cost(world, costs, path)))
x, y = zip(*path)

fig, ax = setup_plot()
fig.set_size_inches(10, 10)
setup_map(ax, world)
ax.plot(x, y)
plt.show()









    



Total path cost: 197

Diagonals and best-first search are considerably worse because they ping-pong back and forth. Once again, this could be the reconstruction algorithm at fault.



In [38]:

    
path = greedy_search(world, costs, reachable, (0, 0), (26, 26))
print("Total path cost: {0}".format(get_total_path_cost(world, costs, path)))
x, y = zip(*path)

fig, ax = setup_plot()
fig.set_size_inches(10, 10)
setup_map(ax, world)
ax.plot(x, y)
plt.show()









    



Total path cost: 128

Slightly better than lateral movements!